home *** CD-ROM | disk | FTP | other *** search
/ BBS in a Box 5 / BBS in a Box -Volume V (BBS in a Box) (April 1992).iso / Files / Bus / M / mathsubs.cpt / FFT.C next >
Encoding:
C/C++ Source or Header  |  1988-01-17  |  6.4 KB  |  189 lines  |  [TEXT/pZIP]
  1. #define RBNEXT  0
  2.  
  3. /*      fft2(a, b, m, dx, inverse)
  4.  *      is an "in place" Fast Fourier Transform for input length a power
  5.  *      of 2.  It will forward transform a time series of real values to
  6.  *      a spectrum of complex values or inverse transform a complex
  7.  *      spectrum to a time series of complex values.
  8.  *    Written by Dale Carstensen, Antares project, Los Alamos National
  9.  *    Laboratory, March 16, 1981 for Unix version 6.
  10.  *    =op operators changed to op= operators for newer than version 6
  11.  *    C compilers by Dale Carstensen, May 14, 1984.
  12.  */
  13.  
  14. double  pi;
  15.  
  16. fft2(a, b, m, dx, inverse)
  17. float   a[];            /* real values of input and returned output     */
  18. float   b[];            /* imaginary values of input (fft2 will zero for
  19.                          * forward transform) and returned output.
  20.                          */
  21. int    m;        /* 1 subtracted from the log (base 2) of the length
  22.              * e. g., 7 for length 256, 9 for length 1024.
  23.              */
  24. float   dx;             /* sampling period              */
  25. int     inverse;        /* zero for forward transform   */
  26. {
  27. int     i;              /* index for permutation and, if inverse, scaling */
  28. int     j;              /* reverse binary index for permutation */
  29. int     n;              /* length of operand    */
  30. int     span;           /* butterfly span distance      */
  31. float   factor;         /* for scaling transforms       */
  32. float   temp;           /* for permutation interchange  */
  33. extern  double  atan2();
  34.  
  35. pi = atan2(0., -1.);
  36. n = 2;                  /* get n = 2**(m+1)     */
  37.     while (m--)
  38.         n <<= 1;
  39.  
  40. /*      permute input to reverse binary order so trig functions will
  41.  *      be required in the normal order.
  42.  */
  43.  
  44. revbin(n);              /* initialize reverse binary counter    */
  45.     for (i = 0;  i < n;  i++)
  46.         {
  47.             if ((j = revbin(RBNEXT)) > i)
  48.                 {
  49.                 temp = a[i];
  50.                 a[i] = a[j];
  51.                 a[j] = temp;
  52.                 };
  53.             if (!inverse)
  54.                 b[i] = 0;
  55.             else if (j > i)
  56.                 {
  57.                 temp = b[i];
  58.                 b[i] = b[j];
  59.                 b[j] = temp;
  60.                 };
  61.         };
  62.  
  63. /*      Perform butterfly calculations using a power of 2 decimation in
  64.  *      time algorithm.
  65.  *      Reference:      G-AE Subcommittee:  The Fast Fourier Transform
  66.  *                      IEEE Transactions on Acoustics, etc. (then Audio
  67.  *                        and Electroacoustics), vol. AU-15, No. 2, June
  68.  *                        1967, pp. 49-50.
  69.  */
  70.  
  71.     for (span = 1;  span < n;  span <<= 1)
  72.         stepfft(a, b, n, span, inverse);
  73.     if (inverse)                /* scale by 1/n */
  74.         factor = 1. / (n * dx); /* scale by 1/n */
  75.     else
  76.         factor = dx;
  77.     for (i = 0;  i < n;  i++)
  78.         {
  79.         a[i] *= factor;
  80.         b[i] *= factor;         /* should be zero, anyway ??    */
  81.         };
  82. }       /* end of fft2(a, b, m, dx, inverse)    */
  83.  
  84. /*      revbin(initialize)
  85.  *      provides a counter in reverse binary order for 0<=count<initialize.
  86.  *      One call with initialize identifies the high order bit.  Each
  87.  *      subsequent call with initialize = 0 (RBNEXT is defined to be 0
  88.  *      for this purpose) will increment the count by one and return it
  89.  *      as an integer function value.
  90.  */
  91.  
  92. revbin(initialize)
  93. int     initialize;
  94. {
  95. static  int     counter, top_order;
  96. register        int     bit, c;
  97.  
  98.     if (initialize)
  99.         {
  100.         top_order = initialize >> 1;
  101.         counter = initialize - 1;
  102.         }
  103.     else
  104.         {
  105.         bit = top_order;
  106.         c = counter ^ bit;
  107.             while (!(c & bit) && bit)
  108.                 {
  109.                 bit >>= 1;
  110.                 c ^= bit;
  111.                 };
  112.         return(counter = c);
  113.         };
  114. }       /* end of revbin(initialize)    */
  115.  
  116. /*      stepfft(a, b, n, span, inverse)
  117.  *      performs one step of the decimation in time butterfly algorithm
  118.  *      for length a power of 2.  The input is assumed to be permuted
  119.  *      to reverse binary order, so the cosine and sine factors can
  120.  *      be generated and used in normal order.
  121.  */
  122.  
  123. stepfft(a, b, n, span, inverse)
  124. float   a[], b[];
  125. int     n, span, inverse;
  126. {
  127. float   angle, cosine, dcossin, inccos, incsin, sine;   /* for trig functions */
  128. float   tempi, tempr, termi, termr;             /* for swapping in place        */
  129. int     i, j;                           /* butterfly indices    */
  130. int     twospan;                        /* for loop termination */
  131. extern  double  pi;
  132. extern  double  sin();
  133.  
  134. /*      Cosine and sine functions are approximated by a technique using
  135.  *      second difference relations.  The formulas are:
  136.  *          func((k+1)a) = func(ka) + inc(k+1, func)
  137.  *              where   func = cos | sin,
  138.  *                      inc((k+1)a, func) = dcossin*func(ka) + inc(ka, func),
  139.  *                      dcossin = -4sin(a/2)sin(a/2),
  140.  *                      inc(0, cos) = 2sin(a/2)sin(a/2),
  141.  *                      inc(0, sin) = sin(a),
  142.  *                      cos(0) = 1,
  143.  *                      sin(0) = 0,
  144.  *                      a = pi/span, and
  145.  *                      k = 0, 1, . . . n/2 - 1
  146.  *
  147.  *      Reference:
  148.  *              Richard C. Singleton:  On Computing the Fast Fourier Transform
  149.  *              Communications of the ACM, vol. 10, No. 10, Oct. 1967
  150.  *              pp. 650, 651.
  151.  */
  152.  
  153. angle = pi / span;              /* for last step, this is 2pi/n */
  154. incsin = sin(angle);
  155. inccos = 2. * sin(0.5 * angle) * sin(0.5 * angle);
  156. dcossin = -2. * inccos;
  157. cosine = 1.;
  158. sine = 0.;
  159. twospan = span << 1;
  160.     for (i = 0;  i < span;  i++)
  161.         {
  162.             for (j = i;  j < n;  j += twospan)
  163.                 {
  164.                     if (inverse)
  165.                         {
  166.                         termr = cosine * a[j+span] + sine * b[j+span];
  167.                         termi = -sine * a[j+span] + cosine * b[j+span];
  168.                         }
  169.                     else
  170.                         {
  171.                         termr = cosine * a[j+span] - sine * b[j+span];
  172.                         termi = sine * a[j+span] + cosine * b[j+span];
  173.                         };
  174.                 tempr = a[j] - termr;
  175.                 tempi = b[j] - termi;
  176.                 a[j] += termr;
  177.                 b[j] += termi;
  178.                 a[j+span] = tempr;
  179.                 b[j+span] = tempi;
  180.                 };
  181.         inccos += dcossin * cosine;
  182.         cosine += inccos;
  183.         incsin += dcossin * sine;
  184.         sine += incsin;
  185.         };
  186. }       /* end of stepfft(a, b, n, span, inverse)       */
  187.  
  188.  
  189.